{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports\n",
    "\n",
    "# utility modules\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import re\n",
    "\n",
    "# the usual suspects:\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# specialty modules\n",
    "import h5py\n",
    "import pyproj\n",
    "\n",
    "# run matplotlib in 'widget' mode\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "sys.path.append(os.path.join(os.getcwd(), '../..'))\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from readers.read_HDF5_ATL03 import read_HDF5_ATL03\n",
    "from readers.get_ATL03_x_atc import get_ATL03_x_atc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data='/srv/tutorial-data/land_ice_applications/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_segs(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot a sloping line for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    #define the heights of the segment endpoints.  Leave a row of NaNs so that the endpoints don't get joined\n",
    "    h_ep=np.zeros([3, D6['h_li'][ind].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li'][ind]-D6['dh_fit_dx'][ind]*20\n",
    "    h_ep[1, :]=D6['h_li'][ind]+D6['dh_fit_dx'][ind]*20\n",
    "    # define the x coordinates of the segment endpoints\n",
    "    x_ep=np.zeros([3,D6['h_li'][ind].size])+np.NaN\n",
    "    x_ep[0, :]=D6['x_atc'][ind]-20\n",
    "    x_ep[1, :]=D6['x_atc'][ind]+20\n",
    "\n",
    "    plt.plot(x_ep.T.ravel(), h_ep.T.ravel(), **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_seg_difference(D6):\n",
    "    \"\"\"\n",
    "    seg_difference_filter: Use elevations and slopes to find bad ATL06 segments\n",
    "    \n",
    "    \n",
    "    Inputs: \n",
    "        D6: a granule of ATL06 data, in dictionary format.  Must have entries:\n",
    "            x_atc, h_li, dh_fit_dx\n",
    "        \n",
    "    Returns:\n",
    "        delta_h_seg: the minimum absolute difference between each segment's endpoints and those of its two neighbors\n",
    "    \"\"\"\n",
    "    h_ep=np.zeros([2, D6['h_li'].size])+np.NaN\n",
    "    h_ep[0, :]=D6['h_li']-D6['dh_fit_dx']*20\n",
    "    h_ep[1, :]=D6['h_li']+D6['dh_fit_dx']*20\n",
    "    delta_h_seg=np.zeros_like(D6['h_li'])\n",
    "    delta_h_seg[1:]=np.abs(D6['h_li'][1:]-h_ep[1, :-1])\n",
    "    delta_h_seg[:-1]=np.minimum(delta_h_seg[:-1], np.abs(D6['h_li'][:-1]-h_ep[0, 1:]))\n",
    "    return delta_h_seg\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 13 data files of which 0 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "#ATL06_files=glob.glob(os.path.join(data_root, 'PIG_ATL06', '*.h5'))\n",
    "data_root='/srv/shared/surface_velocity/'\n",
    "ATL06_files=glob.glob(os.path.join(data_root, 'FIS_ATL06_small', '*.h5'))\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "filename=/srv/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20190608191727_10920311_003_01.h5, exception=name 'map_ax' is not defined\n",
      "filename=/srv/shared/surface_velocity/FIS_ATL06_small/processed_ATL06_20191207103705_10920511_003_01.h5, exception=name 'map_ax' is not defined\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "15f9728901b441e68eca02f9877f2586",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"1092\"\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','05']:\n",
    "    for filename in glob.glob(os.path.join(data_root, 'FIS_ATL06_small', f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure();\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_segs(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
